In [1]:

    

import math
import gpxpy
import glob
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as mdate
from matplotlib import rcParams
%matplotlib inline
rcParams['figure.figsize'] = [16.0, 4.0]


    

In [2]:

    

gpx_file = glob.glob('./gpx/*.gpx')[0]


    

In [3]:

    

def gpx_points(gpx_file):
    with open(gpx_file) as f:
        gpx = gpxpy.parse(f.read())
        for track in gpx.tracks:
            print(track.name)
            for pt in track.walk():
                yield (float(pt[0].time.strftime('%s.%f')), 
                       pt[0].latitude, pt[0].longitude, pt[0].elevation)


    

In [4]:

    

# Mean radius of the earth
EARTH_RADIUS = 6371.009

def haversine(lat1, lon1, lat2, lon2):
    """
    Calculate the great circle distance between two points
    on the earth (specified in decimal degrees)
    
    Return (central angle, distance between points in km)
    """
    # convert decimal degrees to radians
    lat1, lon1, lat2, lon2 = [math.radians(x) for x in [lat1, lon1, lat2, lon2]]

    # haversine formula
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    
    a = math.sin(dlat/2)**2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon/2)**2
    central_angle = 2 * math.asin(math.sqrt(a))

    # s = r * theta
    km = EARTH_RADIUS * central_angle
    return (central_angle, km)


    

In [5]:

    

def pace_tracker():
    idx = 0
    pace = None
    ts = [None, None]
    elev = [None, None]
    lat = [None, None]
    lon = [None, None]
    while True:
        ts[idx], lat[idx], lon[idx], elev[idx] = yield (ts[not idx], pace)
        if idx:
            idx = 0
        else:
            idx = 1
        try:
            dt = ts[not idx] - ts[idx]
            central_angle, km = haversine(lat[idx], lon[idx], lat[not idx], lon[not idx])
            if km == 0:
                continue
            pace = dt / km / 60
        except TypeError:
            pass

def calc_pace(pts):
    a = pace_tracker()
    a.send(None)
    while True:
        ts, pace = a.send(next(pts))
        if pace is not None:
            yield ts, pace

def moving_average(a, n=3) :
    ret = np.cumsum(a, dtype=float)
    ret[n:] = ret[n:] - ret[:-n]
    return ret[n - 1:] / n

def plot_pace(ax, epoch, y, title=None):
    x = mdate.epoch2num(epoch)
    ax.plot(x, y, alpha=0.6)
    ax.xaxis.set_major_formatter(mdate.DateFormatter("%H:%M:%S"))
    ax.xaxis.set_major_locator(mdate.MinuteLocator())
    ax.fill_between(x, 0, y, alpha=0.5)
    plt.xlim([np.min(x), np.max(x)])
    plt.ylim([np.min(y), np.max(y)])
    plt.xticks(rotation='vertical')
    ax.set_ylabel('Pace (min/km)')
    if title is not None:
        ax.set_title(title)
    ax.grid()


    

In [6]:

    

def track_pts(track):
    for pt in track.walk():
        yield (float(pt[0].time.strftime('%s.%f')), 
               pt[0].latitude, pt[0].longitude, pt[0].elevation)

def plot_tracks(gpx):
    fig, ax = plt.subplots()
    for track in gpx.tracks:
        data = np.array([pace for pace in calc_pace(track_pts(track))])
        print(np.mean(data[:,1]))
        print(np.std(data[:,1]))
        data = data[data[:,1] < 20, :]
        pace = data[1:,:]
        window = 1
        ma = moving_average(pace[:,1], n=window)
        plot_pace(ax, pace[window-1:,0], ma, title=track.name)

with open(gpx_file) as f:
    gpx = gpxpy.parse(f.read())
    plot_tracks(gpx)


    

    




6.59581428793
1.96315812532






    

In [6]: